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I. INTRODUCTION 


This report describes a system of computer programs for calculating three-dimensional 
transonic flow over wings, including details of the laminar or turbulent flow in the three- 
dimensional viscous boundary layer. The flow field is calculated in two overlapping 
regions: an outer inviscid flow region governed by the transonic potential equation, and 
a thin boundary layer in which the first order, three-dimensional, compressible boundary layer 
equations are assumed to hold and in which the effects of surface heat and mass transfer can 
be computed. A consistent matching of the solutions in the two regions is achieved iteratively, 
thus taking into account viscous-inviscid interaction. The potential flow analysis is carried 
out by the Jameson-Caughey transonic inviscid wing program FLO 27, and the boundary 
layer analysis is performed by a finite difference boundary layer prediction program developed 
at Boeing. For the wing application, interface programs provide for convenient communication 
between the two basic flow analysis programs. The boundary layer program is general in 
nature and can be used in a wide variety of applications in addition to the wing application 
described in the present report. 

The purpose of the present report is to provide a general overview of the system, including 
detailed descriptions of the theoretical analyses embodied in the boundary layer and interface 
programs (the inviscid wing program is described in detail in Reference 1 ). Readers 
requiring detailed information on use of the programs or program structure are referred to the 
Program User’s Document (Reference 2) and the Program Maintenance Document (Reference 
3). 



II. SYMBOLS AND ABBREVIATIONS 


A ) 

coefficients in 

B 

linearized momentum 

c ) 

and energy equations 

Cf 

7"w 

skin friction coefficient - , 0 

1/2 P e Qe 2 

c* 

attachment line parameter = U e ^ / v 

Cp 

P- 

pressure coefficient - , 0 ; 

/^PcjoQoO^ 


aw e 

8z 


also specific heat at 
constant pressure 


f' , f velocity defect variable: f ' = 1 

PeQe 
, pw 

g’ , g velocity defect variable: g = 1 

P e Qe 


h enthalpy = C p T 


h l> h 3 

H 

K 

K 13> k 31 


metric coefficients 


5 * q2 

2-D boundary layer shape factor ; also total enthalpy = h + — 

6 2 


ratio of adjacent intervals 




Ai?j- 1 


curvature coefficients 


fi x , £ z arc lengths in x and z directions 

M Mach number 


p pressure 

Pr Prandtl number 

Pr t turbulent Prandtl number 

q, q velocity magnitude, velocity vector 

q heat flux 

Q e velocity magnitude in outer flow 

r recovery factor (eqn. 6b) 


2 


R5* 

0 m 

5 * 0 

5 * m Reynolds number = 


R 

Qe 5 *m 
v (77) 


s 

arc length along surface 


T 

temperature 


u 

velocity in x direction 


U S 

velocity in outer flow direction 


u e 

outer flow velocity in x direction 


V 

velocity in y direction (normal to surface) 


V x , Vy, V z 

velocity components from potential flow program 

w 

velocity in z direction 


X 

boundary layer spatial coordinate 


Xp, Yp, Zp 

Cartesian coordinates of potential flow data points 

y 

boundary layer spatial coordinate normal to surface 

z 

boundary layer spatial coordinate 


Pi 

flow angle relative to spanwise (x) coordinate line 

Ps 

flow angle relative to outer inviscid flow 


8* 

displacement thickness 



1 f°° I 

3-D boundary layer length scale = — V 

/9u\ 2 

^*m 

(ay) 

5 *s 

OO 

n J (Pe^e P u s ) dy 
J 0 ' ' 


5 *! 

OO 

%eQe { < PeUe ‘ PU) ^ 


5*3 

= PeQe { ^ eWe_PW) ^ 



1/2 


dy 
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7 

K 

V 

v 

v *f\' v *h 

^ef h 

0, 01, 03 
0 h 

0 s ,0 n 
0,1’ 03 


<f> 

P 

a 


r l’ r 3 


r 

X 


ratio of specific heats 

von Karman constant = .41 
normal direction coordinate 
molcular kinematic viscosity 

effective kinematic viscosities in x and z directions 
effective thermal diffusivity 

effective viscosity factors in boundary layer coordinate system 
effective thermal diffusivity factor 

effective viscosity factors in streamline coordinate system 

01 03 

Rd ’ £d 

dummy variable illustrating difference expressions (pg. 15) 
density 

outer layer effective viscisity function 
shear stresses in x and z directions 

/ 9 9 \ Vi 

shear stress magnitude = m + 73 J 
inner layer effective viscosity variable 


Subscripts: 

e edge of boundary layer 

00 far field reference 

ref arbitrary reference condition 

Overbar: 

— denotes quantities nondimensionalized by their values in the far field 

(°°) 

Prime: 


denotes differentiation with respect to 77 



III. VISCOUS-INVISCID INTERACTION PROCEDURE 


The basic sequence of calculations used to obtain matched viscous and inviscid solutions 
consists of an iterative loop in which the inviscid outer flow analysis and the boundary layer 
analysis are performed sequentially, until either convergence (satisfactory matching) is 
achieved, or the maximum number of iterations specified by the user has been performed, 
as shown schematically in figure 1 . In the first iteration, the potential flow is calculated for 
the bare wing shape, and in subsequent iterations, the effect of the boundary layer flow on 
the outer inviscid flow is felt as a modification to the wing shape through the addition of 
the boundary layer displacement thickness 5*. During each boundary layer calculation, 

5* is computed by integration of the continuity equation. After this new 8* distribution 
has been interpolated from the boundary layer computational grid to the potential flow 
geometry input stations, the 5* to be used for the next potential flow calculation is com- 
puted as a weighted average of the new 5* and the 8* used for the previous iteration. Con- 
vergence is recognized, and the iterations are stopped, when the maximum change between 
the new and old 8*, expressed as a fraction of the maximum 8*, is less than the convergence 
tolerance chosen by the user. At this point a final boundary layer calculation can be carried 
out, if desired. This option allows, for example, the use of a fine boundary layer grid for 
the final calculation, while, to save computing time, a course grid was used for calculating 
the interaction. Finally, machine plotting programs can be called on to display intermediate 
or final boundary layer results. 

The procedure described above is carried out by a set of control cards (NOS procedure, file) 
that loads and executes the various programs in the proper sequence and controls the flow 
of data (disc storage files) between the programs. At the termination of the procedure, 
sufficient data files are always saved to allow the iterative analysis to be re-started where it 
left off, if, for example, the user decides the solutions have not yet converged adequately. 
This feature also allows a user to run the procedure one iteration at a time, with an oppor- 
tunity to examine the results of each iteration and to hasten the convergence of the proce- 
dure by choosing a more nearly optimum 5* weighting factor (under-relaxation factor) or 
boundary layer grid configuration for each iteration. 

The detailed convergence behavior of the procedure depends strongly on the internal con- 
vergence behavior of the inviscid analysis program FLO 27 and on the way in which FLO 27 
is used in the procedure. Most importantly, the final convergence of the inviscid solution on 
a fine grid in FLO 27 tends to be monotonic and can be very slow, especially if pronounced 
shocks are present. Partly as a result of this slow convergence behavior in FLO 27, the most 
economical way of reaching a matched viscous-inviscid solution is not to run the inviscid 
analysis to convergence for every iteration of the interaction procedure, but rather to run 
a relatively small number (between 10 and 50) of internal relaxation sweeps in the inviscid 
code for each repetition of the boundary layer analysis (calculations presented in Section VII 
illustrate this approach). Because of the resulting tendency of the inviscid solution not to 
change drastically between boundary layer calculations, the interaction procedure generally 
converges stably with an under-relaxation factor applied to 8* of about 0.5 to 0.8. The 
number of iterations required to achieve a converged solution depends on the nature of the 
flow, e.g., the presence or absence of shocks and/or boundary layer separation, and the 
choice must be made on an individual basis. A case with entirely subcritical flow and no sep- 
aration in early iterations can converge in three or four iterations, while a more difficult 
supercritical case might require eight or ten. 



IV. BOUNDARY LAYER PROGRAM 


The boundary layer program uses a finite difference method to generate numerical solutions 
to the three-dimensional, compressible boundary layer equations in curvilinear, orthogonal 
coordinates for either laminar or turbulent flow. The structure of the program is dictated 
primarily by the parabolic nature of the equations, which allows the solution to be generated 
in a marching sequence, and by the nature of the initial conditions (initial velocity profiles), 
which the three-dimensional equations require along any boundaries where fluid is flowing 
into the computation region. For either the upper or lower surface of a swept wing, this 
would generally mean that initial conditions would be required at boundaries along the 
wing root and along the leading edge. On a wing with a blunt leading edge, however, a simpli- 
fication is available that allows for a direct solution for the leading edge flow, reducing the 
requirement for initial conditions to those along the wing root boundary. On a swept win[ 
with a blunt leading edge that is free of sharp breaks, the outer inviscid flow will contain a 
surface attachment streamline near the leading edge along which an attachment line (line of 
symmetry) boundary layer analysis is applicable (ref. 4). The attachment line analysis 
is analogous to the laminar similarity boundary layer solution at the stagnation point of a 
two-dimensional, unswept airfoil. Because of spanwise flow in the swept wing case, the 
attachment line boundary layer is not necessarily laminar, but can be turbulent for sufficiently 
high values of the parameter C* (ref. 5). A module is included in the program for solving 
the special boundary layer equations applicable to this attachment line flow for either 
laminar or turbulent conditions. 

The numerical method used in the program is implicit with regard to the solution in the coor- 
dinate normal to the surface, and the differencing in the other two coordinates adapts to the 
direction of the local velocity vector in a manner consistent with the zones of dependence 
and influence in the governing equations. The method is general in nature, and can be applied 
in any surface-fitted orthogonal grid for which some mild restrictions on the velocity field 
are satisfied and for which initial conditions sufficient to determine the boundary layer solu- 
tion can be specified. When the attachment line analysis is used along the leading edge of a 
swept wing, initial conditions are generally only required along the wing root, and the working 
boundary layer grid for either the upper or lower surface of the wing must be constructed so 
that the first spanwise coordinate line on the surface lies along the attachment streamline of 
the outer inviscid flow. A boundary layer grid constructed in this way for a typical jet trans- 
port wing is shown in figure 2. The grid shown was produced by the grid generation program 
described in Section V. 

The marching sequence generates the solution for one surface station (or column normal to 
the surface) at a time. In this sequence the marching proceeds from root to tip along succes- 
sive spanwise coordinate lines, or fines of constant z, starting with the attachment line and 
proceeding to the trailing edge. This marching sequence and the types of difference expres- 
sions used for the x and z derivatives lead to the following restrictions on the velocity field: 


w > 0 

(1) 

-hi A x 

(2) 

h3 A z 


6 



at every point in the boundary layer. The program logic is set up to identify a forbidden 
zone in the solution region and to forbid calculation for any vertical column where the solu- 
tion does not satisfy the above restrictions at all points in the layer. Under some conditions 
(to be discussed later) the boundary of such a forbidden zone will correspond approximately 
to a separation line in the flow field. For attached flow over most common swept wing plan- 
forms, using a working boundary layer grid of the type shown in figure 3.1 , the solution is 
not hampered by restrictions imposed by the numerical method. 

For starting the calculations at the wing root boundary, two options are available: 


1. The flow can be assumed three-dimensional. If no fluid is entering the computation 
region through the boundary, the three-dimensional boundary layer equations can be 
solved. Otherwise, initial conditions (initial velocity profiles) must be specified. 

2. An infinite swept wing analysis can be performed along the boundary, in which the 
spanwise (x direction) derivatives are set to zero. 


A rigorously correct treatment would require option 1 with the initial conditions supplied 
by an analysis of the viscous flow in the root region or by experimental measurements. In 
the absence of such definitive information concerning the wing root flow, the only alter- 
native is to start the calculations using a plausible guess for the initial conditions and to 
acknowledge that the solution over some inboard portion of the wing may not reflect the 
real flow. Comparisons presented in ref. 6, however, show that the portion of the 
wing surface influenced by wing root initial conditions is generally small and that the choice 
of these conditions is not usually of great importance. For most wing planforms, it suffices 
to use the infinite swept wing analysis as a convenient wing root initial condition, even 
through it is not a realistic physical model of the wing root flow. 

BASIC EQUATIONS 

The program solves the three-dimensional, compressible boundary layer equations in curvi- 
linear orthogonal coordinates: 


x momentum: 


z momentum: 


pu 

hi 

pu 

hi 


3U PW 3U 3U ry 

— + — — + (pv) — + puwK 13 -pw z K 31 
3x h 3 3z 3y 



pw 

h 3 


3w 3w 0 

— + (pv) — + puwK 31 -pu z K 13 


1 

hi 

J_ 

h3 


3p 

3n 

3x 

3y 

3p 

, 313 

3z 

3y 


(3) 

(4) 


continuity: — 


3_ 

3x 


thermal energy: 


(ph 3 u) + — (ph i w) + hj h 3 


9y 


(pv) = 0 


pu 3H pw 3H 3H 

- + + (pv) 

hi 3x h 3 3z 


3H 3 / 3h\ 3 / \ 

i7 = ^r ef hiu + 37r 1+WT3 ) 


(5) 
(6 a) 
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As an alternative to solving the thermal energy equation in cases with turbulent flow over 
adiabatic surfaces, the program gives the user the option of using an algebraic formula for 
the density: 


^ = 1 + — rM e 2 (l 

p 2 e V Q e 2 / 


In the above equations h] and 113 are the metric coefficients of the surface-fitted coordinate 
system such that arc length along the surface is given by 

ds 2 = (hj dx ) 2 + dz ^ 2 


r = .89 turbulent 
= .84 laminar 


( 6 b) 


and 


Kl3 = 


1 

dhi 

h l h 3 

dz 

1 

dh 3 

hi h 3 

dx 


The thermal energy equation is expressed in total enthalpy form, where 

q 2 

H = CpT + 2 

The boundary conditions applicable to equations (3) through ( 6 ) are 


dT 3T 

u = w = 0; (pv) = (pv) w ; T = T w or — = — 

dy dy 


w 


at y= 0 


u U e ; w ->W e ; T T e as y 


(7) 

( 8 ) 


The stresses r\ and 73 represent the total effective shear stresses, including both the mole- 
cular shear stress and the Reynolds stress, and are expressed in terms of an effective viscosity 
coefficient, as for example: 


du 

T 1 = p "rfi * 

where 

"efi =^1 
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Here v is the molecular kinematic viscosity, which for compressible flow is calculated as a 
function of temperature by the Sutherland viscosity formula: 


_!L - Pref / T \ 3/2 T ref + 110 
^ref P Vref/ T+110 


T in deg. K 


and (f>\, 03 are effective viscosity factors for which three options are available in the program: 

1. Isotropic eddy viscosity model for turbulent flow: 01 = 03 = 0, where 0 is an effective 
viscosity function originally developed for use in two-dimensional flows and adapted as 
described below. 

2. “Non-isotropic” modification for turbulent flow: Separate 0 \ and 03 are resolved such 
that the component of. the flow in the streamwise (outer flow) direction sees 0 as in 1 ), 
but the cross-flow component of flow sees a 0 whose turbulent part is reduced by a con- 
stant factor, generally .4. The expressions used for 0 j and 0 3 are equivalent to: 

0 s ~ 0 

(^n - l)= 0-4 (0-1) 

3. Laminar flow: 0 ] = 03 = 1 - 

The “non-isotropic” modification for turbulent flow is discussed in more detail in ref. 

6 , where comparisons with two experimental flows are presented. Transition may be 
treated in a crude fashion by using the laminar viscosity for some upstream portion of the 
flow and then switching to the turbulent viscosity. As the program has no provision for pre- 
dicting transition, its location must be specified as an input. For turbulent flow, the effective 
viscosity function 0 is an adaptation of the function used successfully by Mellor (ref. 7) 
for two-dimensional flows: 


0 = 0 

^X 2 y 2 

SI) 

in the inner layer, and 

0 = 0 

/ Qe ^ *m 

\ V 

) 

in the outer layer. 


For computational purposes, an alternate form of the above hypothesis is used, based on the 
value of 0 and the resulting shear stress at the previous step in the iterative solution scheme, 
which is described in the next subsection. When the iterative solution converges, the two 
forms of the hypothesis are completely equivalent. The detailed functions used in the alter- 
nate form are: 


y4 

0 = 1+ + 6 93 * n ^ le * nner ^ a y er > 
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where 


and 

Qe 8* 

4> = o in the outer layer, 

p e 

where a is a modification to the usual constant outer eddy viscosity for better predictions at 
low Reynolds numbers: 


a = .2143 + .1378 logjQ R 5* m 


+ .032/logjo R, 


, \2_.00248(log 10 R«* m ) ; 


and 5* m is a three-dimensional outer layer length scale: 


6 * 


m 


/*°° r ■■■> 

1 / / 3u\ 2 / 9w\2 

qJo y [W/ + \9y / 


% 

dy 


The inner layer is defined to end at the value of y where the inner layer 0 crosses and becomes 
greater than the outer layer 0. The resulting composite 0 thus has a discontinuity in slope, 
but the numerical method is not adversely affected. 

The heat flux q is expressed- in terms of an effective thermal diffusivity coefficient: 


where 


9h 


"ef h = W>h . 

and 0^ is related to the effective viscosity 0 through the turbulent Prandlt number Prp, which 
is taken as a constant: 




( 0 - 1 ) 


Pr t = constant, typically 0.9. 

On the attachment line, the conventional plane of symmetry assumptions are made, leading 
to the following set of equations: 
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x momentum: 


(9) 


pu 3u 
hi 3x 



- 1 3p + 3 
hi 3x 3y 


/ 3u 


) 


z differentiated momentum: 


p /3 w\ 2 pu 3 /3w\ 3 /3w\ 

h3 \3z / hi 3x \3z / 3y \3z / 


+ 



) 


-pu 


2 


3Ki3 

3z 



Continuity: 


Thermal Energy: 



pu 3H 
hi 3x 



3 

h lh 3 — (pv) = 0 
“*( 4+UT 0 • 


( 10 ) 

(ID 

( 12 ) 


where the same effective viscosity and thermal diffusivity formulations are used as with the 
three-dimensional equations, and the boundary conditions are: 

(13) 


(14) 

As in the 3-D case, the algebraic density formula (Equation 6b) can be used as an alternative 
to solving the thermal energy equation. 

Infinite Swept Wing Analysis 

When the infinite swept wing analysis is requested, all x derivatives are set to zero, and the 
three-dimensional equations then require initial conditions only at a single upstream station. 

If, in addition, the attachment line analysis is requested at the initial station, no initial 
conditions at all are required. 

Displacement Thickness Calculation 


3w 

9T 

_ 9T 


— = 0 ; 

(pv) = (pv) w ; T = T w or — 



3z 

3y 

3y 

w 



at y = 0 


3w 3W P 

u U e ; — -> ; T -**• T e as y -*■ 00 

e 3z 3z e 


After the solutions for the velocity and density profiles have converged at a given station 
on the surface, an integrated form of the continuity equation is solved for the physical 
displacement thickness 5*: 


9 r ( 

srM 


p e U e 5* -p e Q e 6* j 


) 9 

+ — hi 
I 3z 1 


Pe w e 5 * -P e Qe 5 *3 


)] 


^1^3 Pw v w 


(15) 
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where 


5*3 = 


jk f o °°U u °- pu ) dy 

Jo 


On the attachment line this equation reduces to: 


“ [ h 3 { Pe u e (5* -5*l)] +hi p e 5* - p e U e 3 =hih 3 p w v w ( 16 > 


The detailed difference equations used in solving these equations are given in the Appendix. 
TRANSFORMATION OF NORMAL COORDINATE 

To put the equations in the form used for computation, the y coordinate is transformed to 


y 



where S* m , the outer layer length scale defined previously in connection with the effective 
viscosity formulation, is used as the normalizing length scale, whether the flow is laminar or 
turbulent. This transformation counteracts the effects of boundary layer growth, keeping 
the velocity profiles scaled advantageously in the computational grid at all stations on the 
surface. For turbulent flow, it has proved to be sufficient to terminate the grid and to apply 
the outer flow boundary conditions at 77 = 15. For laminar flow, a smaller value can be used. 

NUMERICAL METHOD 

The basic equations are non-linear and are solved iteratively by successive substitution, in 
which the non-linear coefficients are evaluated from the previous iteration. The sequence of 
calculations performed by this iterative procedure, in which the momentum equations and 
the energy equations are solved alternately, is shown schematically in figure 3. Because of 
the way in which the effective viscosity depends on the previous iterations, at least six iter- 
ations are performed, even if convergence is indicated earlier. Convergence to within a small 
tolerance in six iterations is not unusual, but if the solution is changing rapidly from one sta- 
tion to the next, larger numbers of iterations, often as high as 20 to 30, are needed. Stations 
near separation display the slowest convergence, but even then an iteration limit of 50 is 
usually sufficient for practical calculations. 

Convergence is recognized, and the procedure is stopped when 5* m and the velocity gradients 
du 3w 

— — and — — at the surface cease to change by more than a specified percentage tolerance 
dy dy 

from one iteration to the next. 
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The x and z momentum equations are treated as equations for the velocity components u and 
w, respectively, which are expressed in terms of velocity defect variables: 


f' 

S' 


1 


pu 


PeQe 

pw 


PeQe 


where primes denote differentiation with respect to 17. The continuity equation is used to 
express the (pv) terms in the momentum equations in terms off, f g, and g'. The continu- 
ity equation is thus eliminated and the momentum equations form a set of two third order 
equations for f and g. Since f and g appear only in nonlinear terms involving first or higher 
order 77 derivatives, the linearized equations solved at each iteration are treated as second order 
equations for f ' and g\ 


The energy equation is treated as an equation for the density ratio, or, since the pressure is 
assumed constant through the layer, the temperature ratio 


Pe T 
d = — = — 

P T e 

The numerical method used to solve the linearized momentum and energy equations is 
implicit, with the equations in effect being written only at the grid points on the unknown 
column normal to the surface. The z derivatives are represented by three-point, second order 
upwind difference formulas which take into account unequal intervals in z and which make 
use of known solutions at (i, k - 1) and (i, k - 2) and the unknown column at (i, k) as shown 
in figure 4(a). When the x component of velocity, u, is positive, the x derivatives are repre- 
sented by the same type of 3-point upwind difference formula in terms of known solutions 
at (i - 1, k) and (i - 2, k) and the unknown column at (i, k). For positive u (w is always 
positive) the difference scheme is unconditionally stable. When u is negative, a 2-point, first 
order difference formula in terms of values on the previous spanwise line (k - 1) is used for 
the x derivatives, as shown in figure 4(b). Note that the sign of u is evaluated locally at 
each point in the layer, and that the different x difference formulas for positive and negative 
u can appear at different points along the same vertical column. In terms of its treatment of 
the x derivatives for negative u, this scheme is equivalent to one proposed by Dwyer (ref. 8), 
and by analogy with Dwyer’s stability analysis, the Courant-Friedrichs-Levy condition appli- 
cable for negative u is assumed to be given by equation (2). 

The above combination of x and z difference formulas is believed to be consistent with the 
concept of zones of dependence and influence in three-dimensional boundary layers (ref. 

9), and accurate results have been obtained for test cases (described below) containing both 
positive and negative u. The detailed difference equations are given in the Appendix. 

The change from second order to first order accuracy in the x derivatives whenever u becomes 
negative does not often cause a serious deterioration in the accuracy of the results. For prac- 
tical wing calculations, it is of little consequence, since negative u seldom appears in the type 
of grid shown in figure 2. However, early in the development of the program, an attempt 
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was made to extend the accuracy to second order for negative u by using difference formulas 
equivalent to evaluating 3-point formulas for the x derivatives at (i, k - 1) and (i, k - 2) 
and extrapolating the values of those derivatives to the unknown column (i, k). This scheme 
displayed good accuracy for small negative values of u, but it became unstable long before the 
assumed Courant-Friedrichs-Levy condition (Equation 2) was violated. This proved to be 
an unacceptable limitation in the calculation of flows such as the duct flow experiment of 
Vermeulen (ref. 1 5), discussed later, and the second order x differences for negative u have 
since been discarded in favor of the first order. 

For special situations in the grid, the program uses alternatives to the x and z difference for- 
mulas described above. When proximity to a boundary makes a third point unavailable for 
either the 3-point x difference for positive u (figure 4a) or the 3-point z difference, the 
program uses 2-point, first order differences. When u is positive, and the solution is not 
available at (i - 1, k) (e.g., if an illegal velocity was encountered there), the x difference is 
taken at the last spanwise line on which solutions at i - 1 and i are available, provided the 
condition 


u hj Ax- 
w I13AZ- 


(17) 


is satisfied, where the difference interval Az' goes back to the spanwise line actually used. 
Figure 4(c) shows this difference schematically for the case where the next spanwise line 
back (k - 1) is used. If none of the available difference meshes can be used legally, solution 
is forbidden, and the program proceeds to the next station. 

The result of replacing the x and z derivatives by finite difference can be written for the mo- 
mentum equations: 



and for the energy equation: 


+ A x f" + B x f + C x 


+ A z f" + B z f" + C z 


0 

0 


(IB) 

(19) 


_Pe_\ 7 Qe_ 7 

Poo / Qoo 7 ~ 1 



0h 

d 



+ Agd 5 + B P d + C e 


0 


( 20 ) 


where the coefficients A, B, C in the momentum equations depend on the solution at the 
previous iteration, and the coefficients in the energy equations depend on d, f, and g at the 
previous iteration and f' and g' from the present iteration. 
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The 17 derivatives (primes) of f' and g' and d at an 77 mesh point j are represented by difference 
expressions of the following forms: 


where 


f .» = fj-M-fj-l 

1 Vj + l-Vj-l 


= \ T 7 j + 1 - Jjj / 1 Vrjj-tjj _ 1 ) 




3 


1 

4>j + ‘A = j (*j + 1 + ^j) 
= \ (*j + *j-l) 


The result of writing these expressions at all the 77 mesh points along the unknown column 
except at the inner and outer boundaries is a set of linear algebraic equations. These equa- 
tions, given in detail in the Appendix, combined with the known inner and outer boundary 
conditions on f', g', and d,are solved by tridiagonal matrix inversion. 

In order to resolve velocity profiles near the surface accurately in turbulent flow, while at 
the same time using an economical number of grid points, it is necessary to vary the grid 
spacing, using small spacings near the surface and larger spacings farther out. To minimize 
the effects of unequal intervals, a grid is used in which adjacent intervals are related by a 
constant ratio: 


K = ^i 

Ai7j + 1 


The initial interval A77 j is dictated by requirements of resolving the inner part of the layer, 
and in turbulent flow, the higher the Reynolds number, the smaller it must be. For Reynolds 
numbers typical of full scale wings, A171 should usually be less than .001 and small values of 
K (near 1) thus require large numbers of grid points. Large values of K, on the other hand, 
can degrade the numerical solution through the loss of second order accuracy associated 
with unequal intervals. Thus, for this source of numerical error, just as with any other trun- 
cation error, there is a trade-off between economy and accuracy. Figure 5 shows results for 
a very high Reynolds number flat plate flow, calculated in several different grids with different 
numbers of points. The calculations all started with the same initial conditions at x = 0, and 
the results are compared at the final station at x = 3.05 m. The results with 256 points and 
with 100 points are nearly identical, indicating that these represent essentially converged 
solutions. Even with only 25 points in the layer and K = 1.5, the results are accurate enough 
for most engineering purposes. A good compromise for this flow and for most practical, 
full-scale wing calculations is to use about 40 points and K = 1 .25. 
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To test the behavior of the numerical method for various flow angles with respect to the com- 
putational grid, especially for negative spanwise velocity u, test cases were constructed from 
a two-dimensional flow for which a two-dimensional solution had been obtained previously. 
For these comparisons, the spillway flow measured by Bauer (ref. 10) was used. A 
comparison of the two-dimensional solution and the experimental results is shown in figure 
6, and agreement is seen to be good. The three-dimensional program was used to solve this 
flow in several different slanted coordinated systems (figure 7) in which the flow appeared 
three-dimensional, that is in which all velocity components and their derivatives along the 
two coordinate directions parallel to the surface were non-zero. The known 
two-dimensional solution from figure 6 was used to provide initial conditions along the up- 
stream boundaries of these calculation regions, and the three-dimensional results were com- 
pared with the two-dimensional solution at the downstream boundaries. Differences observed 
were due entirely to the numerical method, since the same turbulence model was used to 
calculate both the two-dimensional and three-dimensional solutions. When the direction of 
marching resulted in a positive “spanwise” velocity u, the differences in 5* and Cf were less 
than 1/2%, as shown in figure 7(a). Marching in the other direction, with negative u, the 
differences depended on the slope angle of the grid, as shown in figure 7(b). Relatively 
small values of negative u (p = 30°) resulted in differences of less than 1%, while at larger 
values (p = 43°) approaching, the Courant-Friedrichs-Levy condition (P = 45°), the difference 
is increased to almost 2%. 

The flow represented by the above test cases is, of course, undirectional, in that the sign of u 
is the same at every point in the layer. In some cases it is important to be able to predict 
flows in which u changes sign within the layer, but it is more difficult to construct defini- 
tive numerical test cases containing such “cross-over” velocity profiles. Only one flow has 
been calculated containing an extensive region of cross-over profiles (Vermeulen’s curved 
duct flow, to be discussed in the next subsection), and indications are that the present 
method yields accurate results in such cases, provided the Courant-Friedrichs-Levy condition 
is met in regions of negative u. 

Since an important application of the present program is the prediction of separation, it 
is important to understand the behavior of the maching-numerical scheme in the neighbor- 
hood of a three-dimensional separation line. The program logic is set up to forbid further 
calculation on any vertical column where it encounters velocity components which violate 
either condition (1) or (2). When solution is forbidden along one column, it is also forbidden 
along any neighboring columns whose difference formulas must reference the forbidden 
column. Which neighboring columns will be affected, of course, depends on which difference 
formulas are applicable, which in turn depends on the velocity field. Thus, a forbidden zone 
can propagate in the solution region in a way which depends on the velocity field in the 
emerging solution and on the alignment and spacings of the computational grid. Whether the 
boundary of such a forbidden zone approximates the location of a separation line in the flow 
field, thus providing a prediction of separation, depends strongly on the particular choice of 
computational grid. An examination of the mechanics of the propagation of a forbidden 
zone into the grid indicates that the present method, and other marching methods like it, 
can predict the location of a separation line only when the grid is chosen in such a way that 
one or the other of the coordinate line families lies at least roughly parallel to the separation 
line to be predicted. 
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An illustration of the above idea is provided by a hypothetical test case, Case A4, used at the 
recent “Trondheim Trials” (ref. 1 1) of three-dimensional boundary layer methods. 

The task in Case A4 was to compute the turbulent boundary layer on a flat plate with a circu- 
lar cylinder of radius a protruding perpendicular to the plate. As initial conditions, flat 
plate turbulent boundary layer profiles with 0/a = .01 were specified along a cross-surface 
4 radii upstream of the cylinder axis, and as a boundary condition, the two-dimensional, 
incompressible potential flow about an infinite circular cylinder was specified. In a region 
extending laterally 3 radii from the plane of symmetry, the participants at the trials were 
asked to compute the flow from the initial plane as far downstream as possible. All of the 
participants chose to compute the flow in Cartesian coordinates, and as a result, the 
boundary of the forbidden zone propagated laterally from the saddle point of separation on 
the plane of symmetiy, as shown by the straight boundary in figure 3.7. This forbidden 
zone boundary approximates the separation line only in the immediate neighborhood of 
the saddle point. In the experimental flow of East and Hoxey (ref. 1 2) which involved 
the same geometry but had a thicker initial boundary layer, the separation line was nearly 
circular from 0°to 90° A more advantageous coordinate system for calculating the flow in 
the neighborhood of such a separation line is a polar coordinate system centered on the 
cylinder axis. When the present method was applied in polar coordinates, the resulting for- 
bidden zone boundary was the curved one shown in figure 8, with a nearly circular segment 
from 0° to 90°. Along this segment of the boundary, the predicted radial component of the 
skin friction was very small, indicating that a separation line in the predicted flow field had 
been approached very closely. 
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V. COMPARISONS WITH BOUNDARY LAYER EXPERIMENTS 


TWO-DIMENSIONAL FLOWS 

It should be regarded as a necessary but not sufficient condition in the validation of any 
three-dimensional method, that whenever the method is used to predict two-dimensional flows, 
the results should be the same as for an equivalent two-dimensional method containing the 
same turbulence model. The present method can be used to calculate two-dimensional flows 
in either of two ways: The program can be run in the infinite swept wing mode with zero 

aw e 

sweep, or the attachment line module can be used with = 0. Both of these modes of 

dz 

operation have been used to calculate a wide variety of the nominally two-dimensional 
flows used as test cases at the 1968 Stanford Conference (ref. 10) on turbulent boun- 
dary layer computation. In all cases, the degree of agreement with experiment was substan- 
tially the same as that displayed at the conference by programs based on effective viscosity 
models similar to the one used in the present analysis. 

In an experiment by Cumpsty and Head (ref. 5), measurements were made in the 
attachment line boundary layer on the leading edge of a swept strut in a wind tunnel. The 
boundary layer was thin enough, and the span wise variations in all flow quantities were mild 
enough that infinite span wing conditions were effectively simulated on the attachment line. 

A range of values of the attachment line parameter C* was produced by setting the strut at 
various sweep angles. At low values of C* the attachment line boundary layer was laminar, 
and at higher values it was turbulent. The attachment line module of the present program, 
running in the infinite swept wing mode, was used to compute the flow in the turbulent 
cases. Agreement between the predicted and measured velocity profiles was quite good, 
the predicted profiles being essentially the same as those predicted by Cebeci (ref. 13) 
using a turbulence model similar to the present one. 

INFINITE SWEPT WING EXPERIMENT OF 
VAN DEN BERG AND ELSENAAR 

For this carefully executed set of measurements (ref. 14), infinite swept wing condi- 
tions were closely approximated on a swept flat plate, with a chord-wise pressure gradient 
being imposed on the plate by a swept wing-like body suspended nearby. The measure- 
ments covered the rear portion of the plate, which sustained an adverse pressure gradient 
similar to what generally occurs on the rearward upper surface of a lifting swept wing. The 
deceleration of the chordwise component of the outer flow velocity is shown by the decrease 
in streamline angle relative to the leading edge (bottom half of figure 9). Within the layer, 
the streamline angles change more rapidly than they do in the outer flow, and when the 
minimum angle occurring in the layer reaches zero, a swept separation line is indicated. Note 
that the magnitude of the skin friction need not go to zero at the separation line. The experi- 
mental flow is seen in figure 9 to separate in this way, and the location of separation is pre- 
dicted reasonably well by the calculations using the simple non-isotropic effective viscosity 
model, while no separation is predicted by the conventional isotropic model. 
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Well upstream of separation, at x = 1020 mm, the non-isotropic model provides good pre- 
dictions of the velocity and direction profiles through the layer, as shown in figure 1 0. The 
predicted direction profiles exhibit a peak near the surface (corresponding to the minimum 
j3] of figure 9), whose existence seems to be confirmed by the experimental profile, 
though the experimental resolution does not show the detailed shape of the peak. Direction 
profiles predicted by the isotropic model do not exhibit a peak and do not agree as well with 
the data. Closer to separation, at x = 1220 mm, neither model predicts the profiles particu- 
larly well. This deterioration of the predictions near separation, however, cannot be attri- 
buted solely to the skewing of the velocity profiles in three-dimensional flows. In the nomi- 
nally two-dimensional, planar flows studied in the past (see ref. 10), the predictions of sim- 
ple effective viscosity models displayed the same sort of deterioration near separation. 

CURVED DUCT EXPERIMENT OF VERMEULEN 

In this experiment (ref. 15) extensive mean velocity measurements were carried out 
on the flat ceiling of a rectangular duct with a 60° bend. The locations of boundary layer 
measuring stations are shown in figure 1 1 , along with patterns of external flow and limiting 
surface streamlines deduced from the measurements. The pressure along the center of the 
duct (streamwise row C of measuring stations) was roughly constant, but a strong radial pres- 
sure gradient resulted in large cross-flow angles in the boundary layer, as can be seen in the 
limiting surface streamline pattern. 

The flow has been computed by the present method using a computational grid aligned with 
the rows of measuring stations as shown in figure 12. Initial conditions from the experimental 
data were applied along the upstream boundary, row 1, and along the outer side boundary, 
row A. Since the coordinate system is rectangular between rows 1 and 4 and polar between 
rows 4 and 16, the curvature of the system is discontinuous at row 4. The solution region was 
therefore divided at row 4, and the solution was obtained in two steps. First, the solution 
was obtained for the straight section between rows 1 and 4, and the resulting velocity profiles 
along row 4 were saved. These profiles were then used as initial conditions for the solution 
in the curved section from row 4 to 16. Along the inner side boundary, row E, the radial 
velocity in the outer part of the layer is outward from row 8 onwards, and as a result, a 
forbidden zone propagated into the grid from that point and occupied just over half 
the width of the flow by row 16. Because of this, no attempt was made to continue the 
calculations for rows 17 through 19. It is possible that an even finer radial grid spacing bet- 
ween rows C and E would reduce the extent of the forbidden zone, but other grids have not 
yet been tried. 

Calculations were made with both the conventional isotropic effective viscosity and the simple 
non-isotropic model described in Section III. Experimental effective viscosity profiles plotted 
by Vermeulen had shown a reduced cross-flow effective viscosity related to the streamwise 
effective viscosity by a roughly constant ratio of about 0.4 (these data, in fact, were the rea- 
son the simple non-isotropic model was programmed in the present method). It was, there- 
fore, expected that the non-isotropic model would result in better predictions than the iso- 
tropic one. The differences between the predictions turned out to be small, however, and if 
anything, the isotropic model shows slightly better agreement with the experiment. The 
predicted limiting surface streamline directions in figure 1 1 show generally good agreement, 
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as do the predicted boundary layer quantities along row C, shown in figure 13. For the 
velocity and direction profiles at station C-I4, shown in figure 14, the predictions of the 
isotropic model show somewhat better agreement with the data then do the predictions of 
the non-isotropic model. Again, as in the van den Berg and Elsenaar experiment, the non- 
isotropic model predicts a peak in the direction profile near the surface, but in this case, 
there appears to be no evidence for it in the data. Thus, in contrast with the results for the 
van den Berg and Elsenaar test case, the isotropic viscosity model provides better agreement 
with the experiment data in this case. 

SWEPT WING EXPERIMENT OF BREBNER AND WYATT 

In this experiment (ref. 16), velocity magnitude and direction measurements were 
made in the boundary layer of a 45° swept wing of constant chord and AR = 5 at several 
angles of attack. The upper portion of figure 1 5 shows the locations on the planform of the 
profile measurement stations. The potential flow at several angles of attack was analyzed 
using one of Boeing’s panel-type influence coefficient potential flow methods. The boundary 
layer grid generation program (described in Section V) and the boundary layer program were 
run for the highest angle of attack for which boundary layer profiles were measured, a = 6.3, 
and the computed results were found to be in generally good agreement with experiment. 

Velocity profile comparisons are shown for inboard and outboard stations at 80% chord in 
figure 16. The small disagreement at the outer edge of the layer represents the discrepency 
between the potential flow analysis and the measured outer flow. (Note that velocity magni- 
tudes are normalized by the far-field value and thus, do not go to 1 .0.) Within the layer, agree- 
ment is good, considering the limitations of the eddy viscosity model at such low Reynolds 
numbers and the fact that the location of the transition front was not reported precisely for 
the experiment and had to be assumed for the analysis. The small reversal in the slope of the 
direction profiles very near the surface is a result of the simple non-isotropic eddy viscosity 
model, which is discussed in connection with the calculations presented earlier in this section 
for the infinite swept wing (van der Berg and Elsenaar) and the curved duct (Vermeulen). 
Comparisons with the data nearer the trailing edge (x/c=0.99) were not possible because the 
analysis predicted separation at about 96% chord (slightly sooner outboard; slightly later 
inboard). If viscous-inviscid interaction were taken into account in the analysis, the predicted 
separation would move aft or possibly not occur at all. The predicted displacement effect 
cannot be compared directly with the experiment because the experimental data were too 
sparse to allow integration for the experimental three-dimensional displacement thickness. 
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VI. INTERFACE PROGRAMS 


In the viscous-inviscid interaction procedure, the communication of data between the inviscid 
analysis program (FLO 27) and the boundary layer program depends on two basic interface 
programs (see figure 1). The first interface, described in detail below, transforms and inter- 
polates the inviscid outer flow velocity vectors as required for input to the boundary layer 
program and, in the process, generates the curvilinear, orthogonal coordinate grid used in 
the boundary layer analysis. The second interface interpolates the displacement thickness 
calculated by the boundary layer program back to the potential flow geometry input data 
points and is also described in detail below. 

WING BOUNDARY LAYER GRID GENERATOR 

The function of this interface program is to convert potential flow velocity component and 
wing geometry data from FLO 27 ( figure 1 7) into a form usable by the boundary layer 
program and to generate the curvilinear, orthogonal boundary layer grid ( figure 1 8 ). The 
program requires that the wing geometry (x p , y p , z p ) and potential flow velocity (V x , Vy, 

V z ) data be provided at data points arranged along rows or rib cuts roughly aligned in the 
flight direction as shown in figure 17 (data points on a rib cut need not be at strictly con- 
stant y p ). To allow for a correct analysis of the leading edge boundary layer in the boundary 
layer program, the boundary layer coordinate system is constructed such that the spanwise 
coordinate line that divides the wing surface into upper and lower surfaces for purposes of 
the boundary layer analysis coincides with the attachment streamline of the potential flow 
solution. This is accomplished by an iterative procedure in which estimated locations of the 
intersections of the attachment line with the rib cuts are refined repeatedly. At each iteration, 
the directions that spanwise coordinate lines in the boundary layer grid would assume if they 
were to pass through potential flow data points are calculated, using the same quintic curve 
fitting routines that are used later in the construction of the final boundary layer grid. These 
spanwise coordinate line directions are then used to convert the potential flow velocity com- 
ponents V x , Vy, V z at each potential flow data point into a velocity component W p perpen- 
dicular to the spanwise lines and a component U p parallel to the spanwise lines. These com- 
ponents are analogous to the velocity components W e and U e being sought for the final 
boundary layer system. By interpolation, the location is found on each potential flow rib 
cut where W p = 0, and this becomes the assumed location of the attachment line for the 
next iteration. 

The new attachment line intersections result in new spanwise coordinate line directions 
which in turn, result in new velocity components W p and U p , and so on. Usually, after 
about five iterations, the procedure converges to within a very close tolerance such that 
W p = 0 (interpolated) at the assumed attachment line location. This is a convergence criterion 
that could, in principle, be satisfied by any streamline on the surface, as the velocity component 
perpendicular to any streamline is zero. The fact that the procedure converges to the attach- 
ment line in particular results from the special nature of the attachment line as an assymptote 
for all streamlines emanating from the leading edge region. The procedure can be satisfied 
only if it identifies a single streamline that intersects all of the potential flow rib cuts, and 
only streamlines in extremely close proximity to the attachment streamline can satisfy this 
requirement. In practice, the convergence is sufficiently strong that an initial placement of 
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the assumed attachment line on the bisector of the rib cut arc length from the upper surface 
trailing edge to the lower surface trailing edge is sufficient to start the procedure. 

After convergence of the attachment line location, the velocity components Wp and U p 
(known now for the potential flow data points) are interpolated chordwise along the potential 
flow rib cuts to the locations where spanwise coordinate lines of the final boundary layer grid 
will intersect the rib cuts. These intersections form an intermediate potential flow-boundary 
layer grid (or PB grid, as shown in figure 19), and the interpolated velocity components are 
referred to as Wpg and Upg. 

The components Wp and Up are much better suited to chordwise interpolation than are the 
original potential flow components (V x , V y , V z ) which have large second derivatives around 
the leading edge. All chordwise interpolations of velocities and geometry coordinates are 
accomplished using a second order interpolation scheme in which the two adjacent known 
points are matched exactly, and two additional known points are matched in a least squares 
sense. 

The chordwise boundary layer coordinate lines are then constructed step by step as straight 
line segments that intersect adjacent spanwise coordinate lines in such a way as to form equal 
angles with spanwise coordinate line tangent vectors (figure 20), resulting in a very nearly 
orthogonal grid as shown in figure 2. During this process of geometric construction of the 
grid, the spanwise coordinate line segments between PB grid points are treated as fully three- 
dimensional curves in space and are represented by explicit quintics x p (y p ) and z p (y p ), with 
the spanwise Cartesian coordinate y p (figure 17) as independent variable. The quintic curve- 
fitting routines used for this purpose produce a curve fit with continuous second derivatives 
across the defining data points (PB grid points), thus ensuring continuity of the curvature 
coefficient K13. The quintic curve-fitting algorithm was also designed to keep curvature 
as small as possible in regions where the defining data favor small curvature. Wings consisting 
of straight tapered segments with nearly straight spanwise generators, connected to adjacent 
segments of localized spanwise curvature, are particularly well represented by the quintic 
spanwise curves, whose curvature remains localized in a manner consistent with the defining 
data. This is illustrated in figure 21, which shows the curvature of a typical spanwise coor- 
dinate line on the NASA F8 wing. 

The metric coefficients hj and I13 and the curvature coefficients K]3 and K3] are derived from 
the geometry of the final boundary layer grid. Second order finite difference formulas are used 
to compute h] and I13 as derivatives of arc length £ x and 2 Z along the x and z coordinate direc- 
tions, respectively: 


hi 


h 3 


3*x 

3x 

3z 


The curvature coefficient K13, the surface-normal component of the curvature of the span- 
wise coordinate lines, is computed directly from the coefficients of the quintic spanwise 
curves. The chordwise curvature K3 \ is computed by first order finite differences as the 
chordwise derivative of the slope, in the local surface tangent plane, of the spanwise coor- 
dinate lines, also determined from the coefficients of the quintics. 
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The final potential flow velocity components W e , U e in the boundary layer system are computed 
in several steps. First, the components Wpg, Upg, at the PB grid points are converted back 
to Cartesian components V x , Vy, V z , analogous to the original potential flow data velocity 
components. These Cartesian components are then interpolated spanwise from the PB 
grid points to the final boundary layer grid points, using the quintic curve-fitting routines 
for interpolation, with yp as the independent variable. The Cartesian components, known 
now at the final boundary layer grid points, are then converted to a component W e per- 
pendicular to the spanwise coordinate line and a component U e parallel to the spanwise 
coordinate line, using the quintic coefficients to determine the coordinate line directions. 


DISPLACEMENT THICKNESS INTERPOLATOR 

To interpolate the displacement thickness 5* from the boundary layer grid back to the poten- 
tial flow input geometry data points, this interface program requires three basic bodies of 
data: 

1 . A boundary layer solution file containing boundary layer solution quantities, including 
5*, at boundary layer grid points. 

2. A boundary layer grid geometry file containing the Cartesian coordinates of all surface 
points in the boundary layer grids, both upper and lower surface. 

3. A potential flow input geometry file containing Cartesian coordinates of the wing geo- 
metry defining data points. 

Before the interpolation can begin, 5* must be known at all points in the boundary layer 
grid. If flow separation was predicted by the boundary layer program on either the upper 
or lower surface, the 5* data (in 1., above) will be incomplete and must be augmented. This 
is done by extrapolation along chordwise boundary layer coordinate lines from the last valid 
6* data value on each line. The user has the choice of linear or constant-value extrapolation. 
The linear extrapolation option is usually the appropriate choice, but in some special cases, 
such as lower surface separation in the cove region of an aft-loaded airfoil (see F8 wing cal- 
culations, Section VII) constant-value extrapolation is preferable. The required output is 5* 
corresponding to each of the geometry input points in 3., above. For each of these points 
the program uses an iterative search and interpolation algorithm to determine a pair of coor- 
dinates x, z in the curvilinear, orthogonal boundary layer coordinate system corresponding 
to the Cartesian coordinates of the point in question (and thus, at the same time, determining 
whether the point is on the upper or lower surface). A standard two-dimensional interpolation 
algorithm is then used to determine 5* at the potential flow data point, using the boundary 
layer coordinates x, z as independent variables. 
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VII. VISCOUS-INVISCID INTERACTION RESULTS 


WING GEOMETRY AND FLOW CONDITION 


To illustrate the capabilities of the viscous-inviscid interaction system, calculations were 
carried out for the NASA F8 research wing at Moo = 0.8 and a = 2°. The planform and 
defining section locations used are shown in figure 22. Because the planform leading and 
trailing edges display substantial curvature, this wing constitutes a more demanding test 
of the grid generation and boundary layer solution numerics than would a more typical 
transport wing design. Although this example is not at all typical of current LFC wing 
designs, the flow was calculated both with and without distributed surface suction to demon- 
strate the capability of the boundary layer program to treat flows with suction. 

The wing geometry used here is the same one for which Nash and Scruggs (ref. 1 7) 
have previously calculated the three-dimensional boundary layer development. Direct 
comparisons with their results, however, are not possible because experimental pressure 
distributions rather than inviscid flow calculations were used in their calculations, and our 
calculations were made for different, off-design, flow conditions. 

The choice of off-design flow conditions was imposed primarily by limitations inherent in 
the inviscid flow program (FLO 27). Possibly as a result of the large sweep angle of the inboard 
leading edge glove, FLO 27 could not be made to converge at the design Mach number of 0.98. 
A lower Mach number of 0.80 was therefore chosen, where it was possible to achieve a con- 
verged solution. However, as would be expected at a Mach number below the design point, 
a low angle of attack of 2° had to be chosen in order to avoid a strong leading edge suction 
peak. At a = 4.68°, for example, local peak Mach numbers above 1 .8 were present, and 
the strength of the resulting shock clearly violated the assumptions inherent in the use of the 
transonic potential equation. 

GENERAL DESCRIPTION OF THE CALCULATIONS 


The calculations were carried out for five complete cycles without suction and four complete 
cycles with suction, in the general iterative sequence illustrated in figure 1, with numbers of in- 
ternal inviscid sweeps in FLO 27 and 5 * under-relaxation factors as shown in table 1 . In both 
cases, a reasonable convergence of the 5* distribution was achieved. Except for slight move- 
ments of the attachment line location with changes in the outer inviscid flow as the system 
converged, the upper and lower surface boundary layer grids remained nearly the same for 
all of the boundary layer calculations and are shown in plan view in figure 23. In all cases, 
the flow was assumed to be laminar ahead of the arbitrarily imposed transition lines shown 
in figure 23. In lieu of suitable initial conditions along the wing root boundary, the infinite 
swept wing option in the boundary layer program was used to start the calculations. In the 
case with suction, the suction distribution was: 


\ w = -5. for 0 < x/c < .5 
X w = -2.5 for .5 < x/c < 1 .0 


where A w 
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The convergence of the interaction in the case without suction is illustrated in figure 24, 
which shows the 5* distribution, as predicted by each of the five boundary layer calcula- 
tions, along one chordwise lower surface boundary layer coordinate line near mid-semi- 
span. A more global view of the convergence behavior can be seen in figure 25, where the 
maximum corrections to 5*, made by the displacement thickness interpolator program 
(Section VI) for each potential flow geometry defining section, have been plotted against 
defining section number, from root to tip (figure 22 shows the locations of these sections 
in plan view). From figure 24, it is clear that there was some oscillation in the convergence 
of 5*, and it appears that convergence would probably have been faster with a lower under- 
relaxation factor, say 0.5. In view of these results, a factor of 0.5 was used after the second 
iteration in the case with suction, but in that case the boundary layer was much thinner, the 
interaction was weaker, and a fair comparison cannot be made. 


Table l.—lnviscid Relaxation Sweeps and 
5 * Under-Relaxation Factors 


Without 

suction 


With 

suction 



FLO 27 relaxation sweeps 

5* 

Under- 

relaxation 

factor 

Iteration 

number 

Coarse 

mesh 

Medium 

mesh 

Fine 

mesh 

1 

100 

100 

10 

0.8 

2 

100 

100 

10 

0.8 

3 

0 

0 

20 

0.8 

4 

0 

0 

30 

0.8 

5 

0 

0 

40 

- 

1 

100 

100 

10 

0.8 

2 

100 

100 

20 

0.5 

3 

0 

0 

30 

0.5 

4 

0 

0 

40 

- 


In the case without suction the optional algebraic (constant total enthalpy) formulas for 
density through the boundary layer were used rather than solving the thermal energy equation. 
To demonstrate that the viscous-inviscid interaction results would not have been significantly 
different if the energy equation had been used, the upper surface boundary layer for the 
final (fifth) iteration was re-calculated using the energy equation. Figure 26 compares the 
5* distributions predicted by the two calculations, and it can be seen that the types of den- 
sity calculation used had very little effect on 5* in this case. In flows where the viscous- 
inviscid interaction is primarily the result of a turbulent boundary layer, and the surface is 
adiabatic, the algebraic density formula is a good approximation and can be expected to 
give results nearly identical to the full energy equation at some saving in computing time. 

For laminar flow, however, the effective Prandtl number is farther from unity, and, in 
cases where it is important to predict the details of a laminar boundary layer, the energy equa- 
tion should be used. 
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DISCUSSION OF RESULTS 


The calculated pressure distributions display the usual features associated with modem, aft- 
loaded, transonic airfoil sections when operated below the design Mach number and at a 
low lift coefficient. As can be seen in figure 27, which shows the Cp distribution just out- 
board of mid-semi-span, the aft lower surface displays a pressure maximum at about 90% 
chord as a result of the substantial aft camber, and, even at this low lift coefficient, there is 
a pronounced upper surface suction peak near the leading edge corresponding to a small 
region of supersonic flow. The curve labelled “inviscid” in figure 27 was calculated for the 
bare wing shape (40 fine mesh inviscid sweeps in FLO 27 beyond the original 10 sweeps 
shown for iteration 1 in table 1), and in comparison with the other curves, it shows clearly 
the effects of viscous-inviscid interaction. In the case without suction, the boundary layer 
displacement effect produced considerable reductions in lift and pitching moment. In the 
case with suction, the effects are in the same direction but are much weaker. 

The final, converged pressure distributions and details of the corresponding boundary layer 
solutions for the upper and lower surfaces, both with and without suction, are shown in 
figures 28 through 31 . For each surface and flow condition, the outer inviscid solution is 
depicted at its inner boundary (i.e., at the displacement surface) in terms of constant Cp con- 
tours and corresponding inviscid flow streamline patterns (these two plots constituting part a 
in each figure). The corresponding boundary layer solutions are depicted in terms of surface 
shear stress directions, contours of constant displacement thickness 5*, and contours of 
constant 5* s (these three plots constituting part b of each figure). The surface shear stress 
direction field in part b of each figure is shown in the form of surface streamlines, i.e., 
curves constructed to be parallel everywhere to the surface shear stress direction. These sur- 
face streamlines are the computational equivalent of experimental surface oil flow patterns. 
Because 5* in three-dimensional flows is such a volatile quantity, responding strongly to 
convergent or divergent cross-flow patterns, and even taking on negative values, it does not 
provide the familiar indication of local boundary layer thickness that it does in two- 
dimensional flows. For this .reason, the integral thickness 5* s is included to give an indica- 
tion of the thickness of the local streamwise velocity profiles. 


On the outboard wing, the isobars become nearly parallel to constant percent chord lines, 
and in the case without suction, qualitively at least, the boundary layer development follows 
a pattern that would occur on an infinite swept wing subjected to the same pressure distri- 
bution. In regions of strong adverse pressure gradient, S* increases rapidly, and the surface 
streamlines turn outboard. On the upper surface, the most rapid growth in 8 * occurs over 
the last 10% of the chord, while on the lower surface, 5* grows rapidly only over the forward 
portion of the cove. Aft of 90% chord on the lower surface, the pressure gradient is strongly 
favorable, such that 5* decreases, and the surface streamlines turn back inboard. As is typical 
of three-dimensional boundary layers with pressure-driven cross flows, streamline curvature 
near the surface is much stronger than the corresponding curvature of the outer inviscid flow 
streamlines. In the case with suction, the suction is strong enough that large areas of negative 
8 * appear. 

On the aft inboard portion of both surfaces, with and without suction, an intricate, highly 
three-dimensional flow pattern develops. In response to the unsweeping of the isobars in- 
board near the plane of symmetry, jhe flow near the surface diverges strongly, carrying fluid 
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away laterally and leaving a region of negative 5*. This divergence effect can be seen in the 
surface streamline plots for both the upper and lower surfaces and is sufficiently strong to 
produce the decrease in 5 * in spite of adverse streamwise pressure gradients. The local inte- 
gral thickness S* s , however, continues to increase in the streamwise direction. Just outboard 
of this region, the laterally displaced fluid accumulates in a region of converging flow, pro- 
ducing large positive values of 5* but only modestly increased values of 5* s . It is to be expec- 
ted, of course, that the flow pattern described above is influenced to some extent by the 
assumed initial conditions along the wing root boundary (the numerical zone of influence of 
the wing root boundary is somewhat larger than the more strict zone of influence derived 
from the differential equations). And it should be remembered that these initial conditions 
were generated, for the sake of convenience, by the infinite swept wing equations, which do 
not represent a realistic model of the wing root flow. The quantitative details of the predicted 
flow pattern immediately adjacent to the boundary are therefore suspect, but experience has 
shown that the basic qualitative pattern tends to occur, regardless of the wing root initial condi- 
tions. Figure 32 shows 5 * contours predicted for a transport type wing using two widely 
different assumptions for wing root initial conditions (the calculations are described in more 
detail in ref. 6). The contour patterns for both calculations are qualitatively very similar. Pre- 
dicted S* contour patterns such as these will be difficult to compare with experiment, however, 
because directly measurable quantities (local velocity profiles) do not show the strong spatial 
variations displayed by 5 *, and 5* itself would be extremely difficult to derive accurately 
from experimental measurements. In most cases, however, the region in question does not 
constitute a large portion of the wing, and flow field details predicted there have relatively 
little influence on the flow field prediction for the wing as a whole. 

An interesting feature of the present test case is that a massive lower surface separation on the 
outboard wing was predicted when the boundary layer was calculated for the inviscid pressure 
distribution calculated for the bare wing shape. The resulting surface streamline pattern, 
shown in figure 33, contrasts sharply with the corresponding pattern in figure 29, where vis- 
cous-inviscid interaction has reduced the pressure gradient, and no separation occurs. 
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APPENDIX A: DETAILED DIFFERENCE EQUATIONS 
IN THE BOUNDARY LAYER PROGRAM 


This appendix lists the detailed finite difference equations used in solving the following par- 
tial differential equations: 

1. 3-D momentum and continuity equations 

2. Attachment line momentum and continuity equations 

3. 3-D thermal energy equation 

4. Attachment line thermal energy equation 

5. 3-D 5* equation 

6. Attachment line 6* equation 

1. 3-D MOMENTUM AND CONTINUITY EQUATIONS 


These non-linear eauations are linearized by successive substitution, and the resulting linear 
equations are solved in subroutine SOLVEL. The x and z momentum equations are, 


respectively: 

P 

hi uu x 

+■ ~~ wu, 
h 3 z 

+ (PV) Uy 

+ p uwKq3 

- pw2k 31 + ~ (p^efj u y) 

y = o 

© 

© 

© 

© 

© 

© 


£ UW x 

+ ~ ww 7 
h 3 z 

+ (pv) Wy 

+ p uwK 3 1 

-pu 2 K 13 + & -(pn e f 3 w y t 

■< 

II 

O 

© 

© 

© 

© 

© 

© 



where the terms are numbered for future reference in defining the individual coefficients. 
The equations are treated as equations for f' and g' respectively, where: 


f ' = 1 - 




pu 

PeQe 

pw 

PeQe 


and are then expressed in the following computational form: 

“ p Q 2 Kfi u y) y + A x f + B x f + c x = 0 

- ~^2( py e f i w y ) y + A z 8" + B z g' + C z = 0 . 
P e Qe ' 1 7 
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where the coefficients A, B, and C contain the results of expressing the x and z derivatives 
as finite differences and are expanded below. Expressing the 77 derivatives (primes and y deri- 
vatives) as finite differences results in a tridiagonal system of equations: 


-Ip [dO-O] j +A x f" + B x f' + C x = 0 
-|p [d(l-g')] j'+A z g" + B z g’ + C z = 0 


a lj f 'j+1 + a 2j f j + a 3j f j-1 + a 4j = 0 
b lj g j+1 + b 2j g'j + b 3j g j-1 + b 4 j = 0 


where the individual coefficients are 

b u 


b 2 j = 


b 3j = 


Vj+\ 


Vj+\ 


J fej+i+fojjdj-H ] 

-’ij-iL Z J 

■1 [ tej+1 + ±3j) + (?3j + ^3j - 1 ) ] d . + B 

-Dj-i LC’ij + 1 - , jj) 1 ) J J ‘ 

fcj + ^3j-l)dj-l , 


17j+l -nj -1 L (’Jj-’Jj-l) 
-1 


b 4i i)j+i -7)j.| 

+ teiJj3H ) 


tej+i + ^3j) A f to+i + j3j ) ^ toj+^3j-i )l 

| (r)j+l-i?j) (’Jj-’Ij-l) j J 


L (’ij+i -vj) dj+1 


d M 


+ c. 


a li 


1 


’felj+l + £lj) d j+l 


j nj +1 - ’Tj-l l(Vj+\ -Vj) 
-1 


+ 


3.7 • 


"to+i + ^ij) + toj + ^ij-i)" 


J hj+l - T?j-i L(r?j +1 - r?i) (Vj - I7j-1) J 


dj + B x 


a 3 j 


1 

r?j+l -Vyl 


~ felj + ^lj-l) d j- l 


-A- 


a 4j 


-1 


r( 0 ij+i + 0 ij ) 

W “^j-i \_(Vj+\ _7? i) 

.,1 


dj+l 


+ (|ir^iH) d . 


(’ij-’jj-O 


+ C, 


(( gij+i + $u ) feij + ^ij-i )h 

((^j+l-^j) (’Jj-’Jj-l) | J 
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and where 



* 03 

& - Rd 


and R = 


Qe^*m 

v(V) 


The coefficients A, B, and C are sums of contributions from the terms in the differential 
equations: 

A x =A xl +A x2 + A x3 

B x = B xl + B xx 1 + b x2 + b xz2 + b x3 + B x4 

C x = C X 1 + C xx j + C x 2 + C xz 2 + C x 3 + C x 4 + C X 5 + C x £ 

A z = A zl + a z2 + a z3 

B z ~ B zl B zxl + b z2 B zz2 + b z3 + B z4 

Cz = C z ] + C zx i + C z 2 + C zz 2 + C z 3 + C z 4 + C z 5 + C z 5 » 

where the number appearing in each subscript designates the term responsible for the contri- 
bution. The (pv) appearing in terms (3) are replaced by the following expression, derived by 
integrating the continuity equation: 

(pv) " hpH3 [ ~ 5 * m [ h 3%Qe^ - f)] x - h 3P e O e 5 *m y ( *7^ “ 0 - 5 *m [ h 1 P e Qe (7?_g) ] z 
- h lP e Q e 5* r n 2 (i?g'-g) + hih3 (pv) w J . 

The individual contributions are: 

A xl = hy 

Bxl - d'5*m x >1 (2 - r> [dQed - 0] x 

®xxl^ +c xxl = *m x V 

c xi=fe&QeO-n] x 

a x2 = **m z V 0 - s') f 
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b x2 = s* mz n(i-s')f 
5 * 

B xz2 f ' + C xz2 = 0-*')[dQe(l-f)] z 

c x2 = 8*m z >?(l -g') 

d i ^*m r R* 

A *3 " hihij ^ ["3FeQ e 0?-0] x + h 3 (rjf-O + ^ 

[hi^Q e (i7-g)] z + h lS*m z (rfg'-g)- (pv) w j 

6x3 = h]h 3 [ h 3PeQe(^-0] x +h 3 5* mx (^f'-f-rp) 

+ [ h « <* ■ - 4 + h 1* v («' - *) - (p-)w 

Cx3 = + hfhi [h 3 p-e<3e Ol- f)], +h 3 6* mx f- 

[h lp -eQe(„-g)] z -hi «*m z (,?g '* 8 ) + ^ ^>w| 
b x 4 = -g') 

c x4 = 8* m K ]3 d(l -g') 

C x5 = -i* m K 3I d(l-g')2 

c -Ilm! D e°e x W e U e __ _,) 

x6 ~Qe 2 \ hf^ - — - K 13 UeWe + K 31 We 2 | 

A zl = Yi S *m x Vd - f) 

B zl=^y S*m x >7 

Ezxlg ' + C «l =EI^ (l-f)[dQ e (l-g')] x 
c zl = lq J^Ki-n 
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Az 2 = h5 d5 V ! T,(1 “S') 

Bz2 | = r 3 d ' 6 *nV>(2-s')-§^ e feo-s)],) 


^zz2S + t -'zz2 hoO C^Qe ( ^ S )] 


^3Qe 


C z2 “ -jJ 3 d ' 5 ' i: m z 7 ? 


^z3 ^x3 


^ z 3 h]h3 
x* 

L ° m 
PeQe 
_ h I h 3 
PeQe 


[h 3 PeQe(»- 0 ] x +h 3 «* mx (r)f'-f) 

[hip^Q e (T?-g)] z + ,l i 5* mz 0?g'-g -r?) 

(pv) w ! 


*~z3 


d' 

h 1 h 3 

° m 
peQe 


(- [ h 3 peQe <’) - 0] x - h3 5*m x (pf - f) 

r "l hilu 

[h|peQe (y -g)J z + 1)| 5*m z g+ (pv) w 


B Z 4 = -«*mK 3 ld(l - f) 

Cz4 = S* m K 31 d(l -f) 

W.W. - , 

- K 31 W e U e + K 13 U e 2 


C Z 5 = -6* m K 13 d(l -f')2 


^z6 


5*m 

Qe 2 


u e w e 


hi 


x _ 


Replacing the x and z derivatives by finite differences results in the following expressions, 
depending on which adjacent points are used in differencing. Here i is the spanwise (x direc- 
tion) index, and k is the z direction index, with i, k being the present station (unknown column). 
In cases where one index or the other is deleted, it is assumed that that index is not being 
incremented. 
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IZDIFF = 1: 


Ok- 1 


5 *m k ' 5 *m k _i 


m. 


D 


zm 


and similarly for [hjp e Q e (i7-g)l z , [clQ e ( 1— g'>3 z > ^ e z > and ^ e z 

f d6*, 


R xz2 

Cxz2 

R zz2 

^zz2 


m 


^3^zm 


0-g') 


m 


h3Qe 


(1-g) 


J k 

( [dQelk-bQed-f'Jk-l 

^ I D 2m 


m 


h3^zm 
5* 


m 


h3Qe 


[dQe] k - [dQe(l-g')]k-l 


D 


zm 


IZDIFF = 2: 


£ k-2 

k- 1 

k 

z 


m, 


( R zmm " R zm) 5 ; "'m k “ R zmm 5 *m k _i + R zm ^ 


and similarly for [hj p e Q e (r7-g)] z , [dQeCI-g')] z » Ue z j and We z 


R xz2 


^h 3 Qe ^ £( R zmm ' R zm)( dQe) k 

c xz2 = [h 3 0e (1_g ) ], { ( Rzmm ' R zm)( d Qe) k 

k ) 
- R zmm [dQ e O 'f'Jk-l +R zm [dQ e ( 1-^)1 k -2 

B zz 2 = [( R zmm “ R zm) ( -d Qe) k ] 


m k-2 
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Czz2 ^3Qe |(^ ZTTlrn * ^zm)C^Qe] ^ ' I 
tdQ e (l-g')] k _, + R zm [dQe(l-g')] k . 2 


IXDIFF = 0 : Infinite swept wing option 


S*m x = «e x = W 6x = [dQ e (l-f')] x = [h 3 P e Q e (,,-f)] x = 0 


B xxl “ ^xxl B zxl ^zxl B 


IXDIFF = 1 


i-1 i 




and similarly for [dQ e (l-f )] x , [h3p e Q e (i?-f)] x > u e x and w e 


B xxl = 


^-xxl 


l >m l P d Qe 1 

LhlQeJi L DxmJj 

[ 8*m1 [ [dQeJ i ~ [dQ e ( 1 - f)] j - i 

Lhl Qe-Ji ' ^xm 
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L ^1 Qe -»i LDxm J i 


B ~> = ut^ (, - f) J 
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[dQji-CdQeO-g'JH 


IXDIFF = 2 : 


i - 2 i-1 i 


S*m Y = Rxmm-Rxm 


) 5 *m j- 


R X* + R X* 

R xmm ° T R xm ° 1 x 4.2 


and similarly for [dQ e (l-f')] x , [h3p e Q e (i?-f)] x > u e x » and W 6x . 


B xxl “ 


[ B xmm “ B xm) L dQ e J i 
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^'xx 1 “ yijQ e J. |(^xmm ~ ^xini^Qelj ^xmm 

[dQ e (l-f)]i-i + R xm [dQ e ( 1 _ f f )] i-2 } 

B zxl = j^h j Qg 0~f )ij j( B xmm ~ ^xm) L — dQ e J i | 

c zxl = jhjQg j( B xmm “ R xm) t d Qeli 

- B xmm [ d Qe 0“6 )] + ^xm [ d Qe(l”& )] i-2 
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Xxl “ LhlQeJi,k 1 D xp J 

B xxl = 0 

B zxl ~ 0 

_ _ fen n J] | [dQ e (l-g , )3 ii k-1 ~ [dQeO-g^li+^k-l 
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m v = 


r * _ r * 

m i,k- 1 d n ij - !, k - 1 


D 


xm 


and similarly for all other x derivative terms. 

In the foregoing difference formulas, difference intervals with the following definitions were 
used: 


^zm - Zk - i 

Zk " i 

Zm ^k " z k - 2) ( z k - 1 ' Z k _ 2 ) 
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2. ATTACHMENT LINE MOMENTUM 
AND CONTINUITY EQUATIONS 

These equations are solved in a form directly analogous to that used with the 3-D equations. 
The x-momentum equation and z-differentiated momentum equations are, respectively: 


§*m , . §* m P 5* m _ 5* 1 

p e U e 2 P " ef “ y Y + Pe Ue 2 hj ™ x + p e U e 2 pV U V + - tt ?. . Px - 0 
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Pe^e 2 ^ 1 

© ® 
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Making the substitutions: 
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PeUe h 3 
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f' = 1 - 


pu 


g Z = - 


Pe^e 
9w 
P ~ 
9z 




Pe 


3z 


we have the following equations: 


p w. 


peW e 


- [dO-n]'}' +A x f" + B x f' +C X = 0 

- (dg z ') '}' +A z g z " + B z g z ' + C z = 0 
Replacing 77 derivatives by finite differences leads to the tridiagonal system: 


a lj f j + 1 + a 2j f 'j + a 3j f j - 1 + a 4j = 0 

b lj Sz j + j + b 2j gz'j + b 3j Sz'j _ 1 + b 4j = 0 , 
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where 


a = 1 [ l ^)j + l + Wjl d j + l + A 1 

j Vj+\-Vj-\ L (r?j + l-r?j) J 

= -i ri Wj+i + Wj i + i Wj+Wj-i ii d +Bx 
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i mmzi\ d; _ , — a z 1 

r ?j + i- ? ?j-iL( T ?j" 7 ?j-i) J 


^4j ~ C z 

Again, as in the 3-D case, the coefficients A, B, and C are sums of contributions from the 
various terms in the differential equations: 


Ax = A X j + A x 2 

B x ~ B X j+B xx ]+B x 2 

C x ~ C x i + C xx q + C x 2 + C x 3 

A z = A z2 + A z3 
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B z - B z i + B zx 2 + B z2 + B z3 + B z4 
C z = C zx2 + C z5 + C z7 

The continuity equation is eliminated by making the following substitution for Gov): 


(PV) = {"®* m th3Pe u e(’)-0] x - h 3 p e U e S * n , x („f' - f) 

+ 5 *m h lPe u e8z + h l h 3(P v )wJ 


The individual contributions are: 


Axl = h^ d 77(1 _f 

1 8* 

Bxl=^ d ' 5 *m x ’l(2-f')-j^- [dU e (l-f')] x 
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C xl = d ' 5 *m x » 
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5 m hjg z ^ ^ Pv wa u 
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h l 
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S* 

B zl - <W 


a z2 = ^ 

b z 2 = 17j **m x T ' (1 ' f ' ) 


b zx2 8’z + c zx2 “ (1 -f) (dU e g 2 '). 


A Z 3 = — — j J~- [li3Pe u e (v ~ f )] x + I130* 

“ h,li 3 |p e U e 


, .... , , P v wall 

~ h l° mSz " h l h 3 PeUe 
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h l h 3 
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B z4 = -K 31 5 * m d (1 " f ) 


r - 6 * m (w ) 2 - fw \ 

C ^5 - ~ h 3 U e - W h!D e l Wz /x 


^*m K 3 i/tt7 


U e 


Kz ) 


+ S* m Kl3. 


C z7 = - 5 * m Ki 3 z d (1 - f ') 2 

Replacing the x derivatives by finite differences results in the following: 


1 = 1: Infinite swept wing option 

[dUed-OJx =0 


(17 f' - 0 
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[h 3 p e u e = 0 

B x xl = ^xxl = b zx 2 = ^ zx 2 = 0 


5 * mx - 0 


K)x = 0 


1 = 2 : 




i= 1 i = 2 


5 * 


m v ~ 


x * _ js * 

0 m2 mi 


D 


xm 


and similarly for [dU e (1 - f')] x > D^Pe^e ^ " 0] x > ant * (^e z ) x 


B x xl 


^2 

2 h> X m 


C X xl 


S*m2 [dU e ]2'[ dU ed-f')]l 


h 1 2 U e 2 


Dxm 


Bzx 2 


^2 * 5*1112 
1 2 ^xm 


( 1-02 


r _ R & U g g ?'] 1 
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-o ® 


i - 2 i-1 i 


) ^*mj " R xmm k*mj_i + R xm 5 *mj _ 2 


6*m x “ ( ^xmm B. 
and similarly for [dU e (l - f')] x > [ h 3PeUe (V " 0] x > and ( W e z ) x 
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B xxl = 
^xxl = 

b zx2 = 
^zx2 = 


[( R xn'm-Rxm)[-dUe]i| 

[hM ] i { ( Rxmm " R * m ) ^ i 

" B xmm [dUe ( 1 * f )] i - 1 + B xm [dUe( 1 - f )] i - 

'[hiC“e (1 * f) ]- |(^ xmm ’^ xm ) [d^elij 

- ( 1 - Oj . { -Rxmm [d U e g z ’] i - 1 

+ Rxm [dU e g z ]j-2| 


In all of the above expressions, the difference intervals have the same definitions as in 
3-D case. 



3. 3-D THERMAL ENERGY EQUATION 


The equation in total enthalpy form is: 


{/>"efh h y}y 
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where 



7-1 Pe 2 


h 


JL £e 

7- 1 p e 


d 


For purposes of computation, it is treated as an equation for the density ratio d. Like the 
momentum equations, it is linearized by successive substitution, the linearized equation 
taking the form: 


where 


7 _ 7 

Pe Qe 7-1 



+ Ad' + Bd + C = 0 , 


^h = 


_^h_ 

Rd 


and where the coefficients A, B, and C contain the results of replacing x and z derivatives 
by finite differences. Replacing the y or r\ derivatives by finite differences results in a tri- 
diagonal system: 


aid;, i+ao.d + ao d + a 4 . = 0 
A J J 1 Z J J J - 1 4 
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where 


a l; 


a 2 i 


a 3; 


1 


r?j+ 1 -T?j- 1 

_ 7 _ 7 

Pe Qe 7 - 1 

^j+l-Pj-l 


1 

*?] + I '*?]•- 1 



a 4j C 

The coefficients A, B, and C are sums of contributions from individual terms in the differ 
ential equation: 


A = A 9 + A 3 + A 4 
B = B 2 + B 3 + B x 3 + B 4 + B z 4 
c = c x3 + c z4 + c 5 +c 6 
The individual contributions are: 


a 2 

b 2 


( Pe* Oe ~~T + 7 M „2 r e Q e 3 [ 0 -f ') 2 + (l -S') 2 ! d 
P e Q e t 7-1 

^ (7Moo 2 PeQe 3 [0 -O + d -g')] d) 

PeQe l J 


(pv) 


= — f- j - Tf 3 [h3P e Qe (»I - 0] x ' h 3$*m x Ojf - 0 

h l h 3 l PeQe 

-=?T [blP^Qe^-s)] z- ll l 6:!: m z ( 7 ?s'-g) + ll l h 3 777 } 
P e Q e p e^e J 


where 



a 3 


d(_ 7 — . itSV(l-t')-rM» 2 r e Q e 3(l-f')<(S* I11 

h l { Pe Qe y - 1 x 

[(l-f ) 2 + (l-g') 2 ]d| 


1 __ 

^x3^ + ^x3 ~ PeQe^*m 


7 

7- 1 



Bo = +— 
3 hj 


7 Moo 2 p-^Qe 2 5 * m (1 - f’) [(1 - f') |dQ e (l - f')j x 

(1 -g’) fdQ e (l-g')J 


x + ri- «*m x Qe !(l-Of +(l-g)g |d 
6 m x 



A4 = h^|p~ e Q e 771 0 ~ 7 Moo% e Qe 3 ^*m z 0 -g') 

[d -f ') 2 + (l -g') 2 ]d| 


B z4 d + C 24 - ± p- e Q e 8* m ^ (1 -s') { 5 - T ' 1 d ) 2 

B 4 -■! I7M00 2 PiQe 2 5 *m U -g’) [(l-f'){dQ e (l-f)j z +a-g'){dQ e 

d-g')} z + ^ S *m z Qe{(l-f')f" + (l-g')g"}d]j 

C 5 = -7M M 2 peQe 3 {(1 -nrj} 

Cg = - 7 Moo 2 peQe 3 | ( 1 - g') 73 J 


Replacing the x and z derivatives by finite differences results in the following expressions: 

IZDIFF = 1: o k - 1 

® k 
z 
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, > |dQ e (l-f')l k -|dQe(l-f')|k-l 

|dQ e (l-f)} z - 5 : 


'zm 


and similarly for |dQ e (l - g')| z and |h ]P e Q e (.1 - g)} z • 

B z 4 = |^Pe Qe ^*m ^ * 8 ^Jk t° e 1 

c z4 = -[peQe5*m 0 -8')] k [^e d ] k _ , 


IZDIFF = 2: 


7 

A k-2 
A k- 1 


- (Rzmm - Rzm) {dQe U - O } k - Rzmm jdQe ( ■ - O 
+ R Z m{dQ e (l-f')} k _ 2 . 


ldQ e (l -f')}. 


and similary for |dQ e (1 - g') J z and {h\p e Q e (v - g) } z . 


k - 1 


l 


B z 4 PeQe^ + m 


7 


7- 1 
7 


C z 4 - j^PeQe^ m 

" Rzm [p e dj 


7- 1 
* — 7 - 1 


(l _ g)J ^R-zmm - R Z m)|pe J 

(1 ~ g)J R-zmm | p e dj 


k 
k- 1 


IXDIFF = 0 : infinite swept wing option 


{dQ e (l-f')j x ={dQ e (l-g’)j x = {h 3 p- e Q e 07 - 0} x = 0 
Bx 3 = *^x3 = B 
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IXDIFF = 1: 


i- 1 i 
-o ® 


x 


(dQ e (l-nl J dQea-nii-idQeO-nli-i 

l Jx D xm 

and similarly for { dQ e ( 1 -g')} x and {113 p~ e Q e (r? - f)} x • 

B x3 = |VeQeS* m ~ [pj' 1 ]. 

Cx3 = -[P‘eQe«*m ~ (l-f')l. (Ve 7 ’^] 


IXDIFF = 2: 



i- 1 i 

-o 


x 


JdQ e (] - f ) J x (R X inm ~ Rxmj {dQ e (1 - oj. - R xm m {dQ e O ~oj . _ ^ 
+ R-xm jdQe 0 ” f )} | __ 2 5 

and similarly for {dQ e (l -g')} x and {113 p~ e Q e (rj - f)} x . 
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dQp(l 


-n) ; 


= ldQed-f')l j.k-i-ldQeU-f 1 )! i+l,k-l 


D 


xp 


and similarly for {dQ e (l -g')} x and {h 3 p e Q e (i7-o} x • 


b x3 = 0 

^x3 ~ L°e Qe 8*m r H - f )1 

L 7-1 J 1 , k 



IXDIFF = 5: 


i - I i 

k- 1 

-o 6- 



dQe (1 ' f ) 


_ j dQe ( 1 ~ f )} i, k - 1 ~ jdQ e (1 - f )j j - 1 , k - 1 
D, 


xm 


and similarly for { dQ e ( 1 - g')} x and {h 3 Q e (V - 0 } 


x • 



49 



4. ATTACHMENT LINE THERMAL ENERGY EQUATION 

The development here is directly analogous to that used with the 3-D equation. The basic 
equation is: 



+ (pv) H y + —p uH x 
y h l 

© 


-M 

© ® 


= o , 


which results in a linearized equation for the density ratio d: 


-p e ^U e - — ^ dj + Ad f + Bd + C — 0, 

in which 

A = A 2 + A 3 
B = B 2 + B 3 + B x 3 
c = c x3 + c 4 

The resulting tridiagonal system is: 


where 


a *j 


a 2; 


a 3 


1 

^j + 1 -*} 
Pe 7 Ue 7 


T- 1 


bj + l-bj -1 

1 


J i?j + 1 - - 1 


dj + a 3j dj . ] 

+ 

II 

O 

( 

f th i + i + £h\ i4 ' 

’Jj + 1 ' Vj ) 

7 © + 1 + 

* /ihj^hj.A 

Vnj + i-ij , 

/ \ r?j - T7j _ 1 / 

[-^eUe y 7 _ , 

4 ‘ 

\ T7j - 77 j _ 1 / 


3 4j =C 
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The individual contributions to A, B, and C are: 


A 2 =~ |^ 7Ue ^T"l + 7Moo 2 p e U e 3(l-f)2 d 


ee 


b 2 = {7Mc«2p~ e U e 3(l-Ord} , 


P*U, 


ee 


where 


(pv) 

Pe u e 


5 * 


^3 j~P~uf t h3 Pe U e^“^]x-h3 5*m x .(^f'*0 + 5* m h 1 g z 


+ hjh 3 


(pv) 


w 


! 


^e U e 

A3 = ^ {-P e 7 ° e ^ t?8* rnx (l-f')- T M 00 2^- e u e 3^5* mx (i-f)3 d J 


B x3^ ^ ^x3 — — p e U e 5 '' 


m 


7- 1 


(1-0 jp^ -1 dj 


B 3 =2 7 M„ 2 P - Ue 2 8 . ffl (l .02 |[dU e (l -f)] x+ |_S*m x U e df' 

C 4 = -TM„2p- e U e 3 {(1 - O r^j* 


Substituting finite differences for the x derivatives yields: 
I = 1 : Infinite swept wing option 


{dU e (l-f')) x = {h 3 p^U e (t)-o} x = 0 


b x 3 Cx3 - 0 


I = 2 : 



dU e (1 -f) 


H 


|dU e 0 -O| 2 - |dU e (1 - f')} 1 


D 


xm 



and similarly for j I13 p e U e (17 - f) | 


Bx 3 = A e U e S* m (1 - f)J. k r 1 J. / °xm 
Pe U e 5* m ^ (1 - nl k -> d J i _ t / D 


^■x3 


I > 2: 


— o — $5 — *■ 


i - 2 i - 1 i 


{dU e (1 - f )) x = (Rxmm.- Rxm) {dU e ( 1 - f')}. - R xmm { dU e ( 1 - f) j . _ 
+ R xm {dU e (l-f)}._ 2 , 
and similarly for {^PeUe^ - 0 } x 


x 3 


7 - 1 


C 


x 3 


P; U e 5 * m-^T (1 

P e U e 5 *m d 

7~ 1 

r 

+ R 


i xm 




(kmm'km) / 

" f )J . j* ^xrr 

J 


_ 7 -1 

Pe 


Ji 


Pe 7 d 


i- 1 


5. 3-D 6* EQUATION 


The equation is: 


h 3 {Pe u e 5 *-PeQ e 5 *l 


hi p- e W e S*-p e Q e S* 3 


)] 


= hih 3 p e Q e 


— ^ P w v w 


PpQ 


e^e 
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When derivatives are replaced by appropriate difference expressions, the equation becomes: 


[A LX S* + D dx ] + [A LZ S* + D DZ ] = hj h 3 p e Q e 

6 C 


or 


-(DdX + Ddz) +h l h 3 P e Q e ^f 

5 * = / — ~ t » 

^A L x + Alz j 

where Alx> ^DX» ^DZ> and Al Z depend on the types of x and z differencing used as 
follows: 

IZDIFE = 1 : 6 k - 1 

$ k 
z 

\ 

At 7 = [h p W 1 

Dzm L 1 e eJ >, k 

°DZ = ~ { [- h,p e Q e 5* 3 ] i ^ - [h](PeWe5*-p e Q e 5* 3 )] jjk . j 


IZDIFE = 2: 


A k-2 
A k - 1 

k 


AlZ zmm " ^ zm) t' 1 lPe^e ] i, k 

^DZ ~ (^zmm “ ^zm )lr^ l 1 3I i, k " ^zmm Ij’i^e^e^* 

-W*3)]i,k- 1 +R zm Ij 1 1 (pe W e^* " i, k - 2 
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IXDIFE = 0: Infinite swept wing option 


a lx = o 

d DX = 0 

IXDIFE = 1: o <8 

i-1 i x 

a LX = fr— - [ h 3PeU e ]i,k 
xm 

^DX ““ pv | [-li3P e Qe^*l]i, k ' _ p e Q e 6* i) i _ j 5 

xm ' 

IXDIFE = 2: — ° o <8^— x 

i-2 i-1 i 

A LX = ( R xmm- R xm)[ h 3^e U e]i,k 

^DX ~ ( R xmm " R xm ) [j^^eQe^* lj i, k “ R xmm [QlJ i-l,k 
+ R xm [Ql] i - 2, k 

where 

[Ql\k = [ h 3(p e U e «*-PeQe 5 *i)] ijk 

IXDIFE = 3: i i + 1 

k - 1 o o 

k j 'x 

z 

A LX = o 

° DX = { [Qi] i, k-i - Witi,k-iJ 
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IXDIFE = 5: i-1 i 



ALX = 0- 

Ddx = ^- {[Qi 3 i, Jc- 1 - [ Q l] i -l.k-ll 

u xm l > 

Cross-Over Velocity Profile: In the special case when the u velocity profile is of the cross- 

over type (u is both positive and negative along the same column) the zone of dependence 
of the 5 * equation is not properly covered by any of the above difference expressions. In 
this case, a zig-zag difference expression is used when the required points are available: 

i-1 i i + 1 

k- 1 
k 

AlX = R x .pm [ h 3Pe U e] i, k 

d DX = R xmp | [Pi ] i + 1 , k - 1 ’[Q^i, k-lj 

+ R xpm{[- h 3P e Qe 5 *l]i, k ’ C Q l] i- 1 , k} » 

where 

[Ql] k = [h 3 (p e Ue 6 * ' W*l)l j, k 

R = ~ D *P 

xpm " D xm (-D xp+ D xm ) 

R _ D xm 

xm P - D X p (-D X p + D xm ) 
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When the point at i - 1 , k is not available, an alternate expression involving three points 
on the previous spanline is used: 


i- 1 

i i + 1 

k - Q 

f- 


— v: 

X 


A LX = 0 

d DX = R xmp [Ql] i + i ( k - 1 + ( R x P m " R xmp)[Ql]i 5 k _ \ 
- R xpm [Qllj. l,k - 1 


6. ATTACHMENT LINE 5* EQUATION 


The equation is: 


[l’3 ku e ( 6 * - 6*l)}] x + l'l [p e W ez 5* -p e U e 6*3 z ] 


, , Tt P w v w 
hih 3 p e U e7 - u - 


When the x derivative is replaced by appropriate difference expressions, the equation 
becomes: 


[ a lx 6 * + d dx] 


[lMp e W ez 8*-hiP e U e 8* 3; 


_1 u - TT P W v W 
h l h 3 p e U e v 


( h I Pe u e 5 " : 3 Z - °DX ) + h i 11 3 % u e XO 7 

5* = t _ — — v 

( A LX +h lP e W e z ) 

where A lx and D DX depend on the particular x-difference expression, as follows: 
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1=1: Infinite swept wing option. In this special case, the equation reduces to: 


I = 2: 


I > 2: 


8 * 



P w v w 
3 z + h 3 p e Q e 


) 




X 


h3P e U e 

a lx = -d: — 


xm 


d dx - 


D 


xm 


{(-h3P e U e 5*l) 2 - [h 3 P e U e (6^«*l)] 



— o -<& 

i- 1 i 


^LX (^xmm ^xm) [j^Pe^e] i 
^DX = (^xmm -R xm) [" ^3^e^e ^ 1 ll [ 

‘ R xmm CQl] i _ i + R xm [Ol], _ 2 
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Figure 7 . — Viscous-inviscicJ Interaction Procedure 





X 




Figure 2.— Curvilinear, Orthogonal Coordinate System for Boundary 
Layer Calculations on a Swept Wing 
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Figure 3.— Iterative Procedure For Boundary Layer Solution of a Single Station 
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a) Positive u 


b) Negative u 


c) Alternate form, positive u 



column 


Figure 4.— Difference Molecules for Positive and Nt 
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Calculation started at x = 0 Ft ; = 0.152 M 



Figure 5.— Effects of the Ratio K of Adjacent 77 Intervals on the 
Calculation of a Flat Plate Boundary Layer at High Re: 
M = 2.80, Re = 8.07 x 10 7 AT 7 . 



2 3 4 5 6 7 8 9 10 11 12 


x Ft 

Figure 6.-2-D Solution for Bauer's 60° Spillway Flow 
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Figure 11.— Comparison of Predicted Limiting Surface Streamline 
Directions With Streamline Pattern Measured By 
Vermeu/en 


Initial conditions 



Figure 1 2.— Computational Grid Used for Curved Duct Flow 
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Figure 13.— Comparison of Predicted and Measured Boundary Layer Quantities Along 
Central Row (C) of Measuring Stations for Data of Vermeu/en 











Figure 17.— Potential Flow Coordinate System 
With (x p , y p , z p j and (V x , V y , V z ) 

Known at Data Points Arranged 
Along Rib Cuts 


Potential flow rib cuts 



Figure 19.— PB Grid Consisting of Intersections 
of Boundary Layer Spanwise Lines 
and Potential Flow Rib Cuts 



Figure 18.— Boundary Layer Coordinate System 
With Ue, We Known at Grid 
Intersections 



Figure 20.— Final Boundary Layer Grid 
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C L 

C M 

Inviscid 

n an? 

-0.217 

No suction-converged 

0.370 

-0.192 

With suction-converged 

0.391 

-0.209 



Spanwise location 
Yp = 1 3.49 in 


Figure 27.— Predicted C Q on a Section Near Mid-Semi-Span 
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Inf span analysis 
at root 

Quasi 2-D analysis 
at root 



0.1 


u.z. \ 
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slant 5 * Predicted for a Transport Type Wing 
r ent Wing Root Initial Conditions. 
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